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Current  seismic  yield  estimates  are  based  on  scaling  laws  which  relate  the  seismic 
source  spectrum  to  the  yield  and  depth  of  burial  (Mueller  and  Murphy,  1971;  Helmberger 
and  Hadley,  1981;  von  Seggern  and  Blandford,  1972).  Since  these  scaling  laws  do  not 
contain  any  physics  relating  to  the  inelastic  source  regime,  they  must  be  calibrated  for  each 
site  using  shots  of  known  yield  to  fit  source  parameters.  However,  there  is  reason  to 
expect  that  the  source  parameters  in  such  scaling  laws  are  especially  sensitive  to  the  pre¬ 
shot  fracture  spectrum  in  the  source  emplacement  rock.  Non-linear  processes  like  crack 
growth  and  amplitude  dependent  attenuation  are  known  to  depend  on  the  density  and  length 
distribution  of  preexisting  cracks.  By  implication,  both  the  size  of  the  elastic  radius  and  the 
pulse  shape  at  that  radius  should  depend  on  the  fracture  spectrum.  The  ultimate  goal  of  the 
research  reported  below  is  to  formulate  a  source  function  with  parameters  which  can  be 
determined  from  measurable  physical  parameters  of  the  source  rock  such  as  velocity, 
density,  porosity,  initial  fracture  spectrum,  and  degree  of  water  saturation.  Such  a  scaling 
law  would  allow  more  accurate  corrections  for  variations  in  the  emplacement  medium 
without  requiring  extensive  calibration  shots. 

There  is  observational  evidence  that  the  source  function  is  sensitive  to  the  fracture 
spectrum  in  the  emplacement  medium.  Rimer  et  al.  (1987)  performed  numerical 
calculations  of  the  particle  velocity  measured  in  the  free  field  for  underground  tests  in 
granite.  Using  laboratory  data  for  rock  strength,  they  were  unable  to  predict  pulse-width 
successfully.  In  order  to  make  the  models  correspond  to  the  field  data,  much  lower 
fracture  strength  than  those  in  the  lab  were  required.  Sammis  (1989)  offered  the 
explanation  that  the  weaker  rheology  was  the  direct  result  of  larger  fractures  in  the  field. 
More  recently,  McEvilly  and  Johnson  (1989)  have  measured  the  seismic  radiation  from  a 
"'~!cs  cf  chemical  explosions  as  a  function  of  depth  in  a  limestone  quarry.  As  illustrated  in 
Fig.  1 ,  the''  observed  a  pronounced  shift  towam  higher  frequencies  as  die  depth  of  buual 
increased.  We  hypothesize  that  this  is  due  to  the  suppression  of  crack  nucleation  and 


growth  with  confining  pressure.  In  a  recent  review  of  source  models  and  scaling  laws, 
Denny  and  Johnson  (1991)  conclude  that  the  cavity  radius,  the  seismic  moment  and  the 
comer  frequency  of  explosions  are  all  dependent  on  the  depth  of  burial.  They  reference 
Larson's  (1984)  work  which  suggests  that  these  effects  may  result  from  the  increase  in 
shear  strength  with  depth.  Indeed,  since  shear  failure  in  the  damage  mechanics  model  is  a 
simple  consequence  of  the  growth  of  fracture  damage  to  a  critical  level  (where  the  stress- 
strain  curve  becomes  unstable)  the  suppression  of  crack  growth  by  confining  pressure  is 
equivalent  to  an  increase  in  shear  strength  with  depth. 

If  the  nucleation  and  growth  of  fractures  in  the  non-linear  source  regime  can 
significantly  broaden  the  source  pulse  and  change  the  elastic  radius,  then  the  process  may 
significantly  reduce  the  radiation  of  seismic  waves  at  the  high-frequency  end  of  the  seismic 
spectrum.  This  is  especially  significant  since  recent  yield  verification  and  source 
discrimination  schemes  increasingly  utilize  higher  frequency  local  phases  as  the  joint 
verification  program  has  allowed  more  local  seismic  monitoring.  The  spectrum  of  initial 
fracture  sizes  at  a  test  site  may  turn  out  to  be  as  important  as  emplacement  rock  type  in 
determining  the  yield  from  seismic  radiation. 

SEISMIC  .SOURffiJEUMOJQNS 

The  uncertainty  in  seismic  yield  determination  can  be  broken  down  into 
uncertainties  associated  with  the  propagation  of  the  seismic  waves  and  uncertainties 
associated  with  the  source  coupling.  The  recent  use  of  high-frequency  local  phases  such  as 
Lg  (see  e.g.  Ringdal,  1990)  has  significantly  reduced  uncertainties  associated  with 
propagation,  but  our  understanding  of  the  effects  of  the  emplacement  medium  on  seismic 
coupling  has  not  significantly  improved  since  the  1971  work  of  Mueller  and  Murphy. 

The  Mueller-Murphy  source  assumes  a  step  function  followed  by  an  exponential 
decay  of  the  pressure  source  which  acts  on  a  spherical  surface  at  the  "elastic  radius"  rci. 

p(t)  =  (p0e  ' at  +  p0c)  H  (t)  (1) 
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This  time  function  was  chosen  because  it  mimics  near-field  velocity  records.  The  elastic 
radius  is  chosen  such  that  it  contains  all  the  non-linear  source  phenomena;  the  medium  is 
assumed  to  be  perfectly  elastic  for  r  >  re). 

The  elastic  radius  is  defined  by  the  assumption  that  the  peak  stress,  pos,  at  the 
elastic  radius  is  a  fixed  fraction  of  the  overburden  pressure  at  the  d^pth  of  burial  h. 

pos=1.5pgh  (2) 

The  peak  shock  pressure  in  the  non-linear  regime  is  assumed  to  be  of  the  form 


p,=^. 


r"  <3) 

They  further  assume  that  n  is  independent  of  the  emplacement  medium  and  m  =  n/3,  so  the 
final  scaling  relation  is 


rel] 

rel2 


(4) 


which,  when  the  explosions  are  in  the  same  medium,  reduces  to 

lehjwn  1/3  (M  1/0 

rel2  (W 2 1  lh]  /  (5) 

P  wave  amplitude  spectra  in  the  far-field  were  used  to  determine  the  empirical  constants:  n 
=  2.4,  Asait/At-r  =  12,  AshaleMi-r  =  5.3,  and  Avt/At-r  =  0.23  (whure  subscripts  t-r  =  tuff- 
rhyolite  and  vt  =  volcanic  tuff).  The  constant  a  in  equation  (1)  is  assumed  to  be  of  the 
form 


a  = 


(6) 


where  c  is  the  p-wave  velocity  and  k  depends  on  the  source  rock:  for  tuff,  k  =  1.5; 
rhyolite,  k  =  2.0;  shale,  k  ~  2.4;  and  salt,  k  =  4.5.  The  broad  pulses  observed  in  granite 
imply  a  very  small  k,  which  suggests  that  k  may  reflect  the  fracture  structure  as  well  as  the 
porosity.  Using  our  damage  mechanics  model,  it  should  be  possible  to  relate  k  to  the  initial 
damage  spectrum. 
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Hence,  the  Mueller-Murphy  source  makes  a  number  of  reasonable  assumptions 
about  the  source,  and  then  uses  calibration  shots  to  fix  the  values  of  the  media  dependent 
parameters  A,  n,  and  k.  One  objective  of  the  work  described  below  is  to  use  recent 
developments  in  the  damage  mechanics  of  brittle  solids  to  make  a  physical  interpretation  of 
the  coupling  parameters  of  the  Mueller-Murphy  source  model  in  terms  of  physical, 
measurable,  properties  of  the  source  emplacement  medium,  thereby  eliminating  the  need  for 
calibration  shots  at  each  new  site,  and  allowing  an  estimate  of  the  variation  in  the  source 
function  as  a  function  of  geological  setting  at  a  given  test  site. 

To  help  focus  the  discussion,  consider  the  schematic  diagram  of  buried  explosive 
source  shown  in  Fig.  2a.  For  our  purposes,  we  simplify  the  source  and  identify  three  non¬ 
linear  regimes  as  indicated  in  Fig.  2b:  the  "hydrodynamic  regime"  in  which  rock  flows,  the 
"damage  regime"  in  which  the  rock  behaves  as  a  solid  but  stresses  are  large  enough  to 
extend  existing  cracks,  and  the  non-linear  attenuation  regime  (not  shown  in  2a)  in  which 
stresses  are  large  enough  to  produce  amplitude  dependent  attenuation  by  motion  on 
preexisting  fractures  but  not  sufficiently  large  to  cause  additional  fracture.  The 
hydrodynamic  radius,  q,,  depends  on  the  equation  of  state  of  the  emplacement  medium  and 
is  the  subject  of  high  pressure  shock  wave  studies.  The  damage  radius,  rd,  is  defined  by 
the  condition  that  the  peak  radial  stress  has  fallen  to  a  level  which  is  just  sufficient  to 
nucleate  fractures  from  initial  flaws  in  the  emplacement  medium. 

The  damage  mechanics  developed  by  Ashby  and  Sammis  (1990)  allows  a 
quantitative  evaluation  of  rd-  Their  equation  for  the  radial  stress  Gr  required  to  initiate  flaws 
when  the  hoop  stress  is  gq  is 

Gr  =  G0  +  C1Ge  (7) 
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where  ci  depends  on  the  coefficient  of  fnction  p  on  the  starter  flaws 

V  i+ii2  +  (i 

ci  =  . - -  - 

V  l+p2  -  p 

The  other  constant  ao  is  defined  as 

a  = _ 3S _ (Kic 

V7^-pW™ 

where  Kic  is  the  critical  stress  intensity  factor  for  tensile  failure  and  a  is  the  half-length  of 
the  largest  initial  flaws  in  the  emplacement  medium. 

It  is  important  to  note  that  the  damage  radius  is  not  simply  a  function  of  rock-type. 

In  fact,  the  parameters  Kic  and  (i  are  almost  independent  of  rock  type.  Rather,  rj  is  most 
sensitive  to  the  size  of  the  largest  flaws  in  the  emplacement  medium.  The  effects  of  ground 
wa'er  saturation  is  to  reduce  the  effective  p  on  pre-existing  cracks  thereby  increasing  r<j. 
Note  that  eqn.  (7)  is  of  the  form  assumed  by  Mueller  and  Murphy  (1971)  (eqn.2  above)  as 
long  as  the  depth  of  burial  is  great  enough  that  cjae  »  o0- 

The  elastic  radius  rci  is  more  difficult  to  define  because  there  is  no  physical  cutoff 
to  the  non-linear  attenuation.  However,  if  amplitude  dependent  attenuation  is  due  to  motion 
on  pre-existing  flaws,  than  rc|  can  be  expected  to  scale  with  flaw-size  in  a  similar  wav  as 
the  damage  radius  since  a  smaller  stress  is  required  to  produce  motion  on  a  larger  fracture, 
in  fact,  if  the  source  emplacement  medium  is  heavily  jointed,  the  elastic  radius  could  be 
very  large  indeed.  The  possible  role  of  joints  is  under  investigation  by  several  groups  (see 
egn.  Ileuz'e  et  al,  1991 ) 
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The  immediate  objective  of  this  project  is  to  determine  the  change  in  pulse  shape 
associated  with  the  extensive  rock  fracture  which  occurs  between  the  hydrodynamic  radius 
and  the  damage  radius  (see  Fig.  2).  Further  modification  of  the  pulse  caused  by  amplitude 
dependent  attenuation  beyond  the  damage  radius  is  being  investigated  by  other  groups. 
Which  process,  if  either,  is  more  important  remains  to  be  determined. 

A  complete  damage  mechanics  suitable  for  incorporation  into  the  numerical  codes 
which  simulate  underground  explosions  has  two  components: 

a)  A  model  for  the  nucleation,  growth  and  interaction  of  tensile  cracks  which 

relates  crack  length  to  the  applied  stress  field,  and 

b)  A  model  for  the  effective  elastic  constants  as  a  function  of  crack  length. 

Part  (a)  has  been  completed  and  verified  by  predicting  the  fracture  nucleation  and  failure 
surfaces  of  wide  range  of  rocks  as  measured  in  triaxial  laboratory  experiments  (Ashby  and 
Sammis,  1990).  We  are  now  focusing  our  attention  on  Part  (b),  which  is  the  main  subject 
of  this  report. 

Before  discussing  the  elastic  constants,  it  is  interesting  to  see  how  the  damage 
mechanics  may  be  incorporated  in  the  numerical  source  codes. 

1 )  The  current  stress-state  is  used  in  the  equations  of  motion  to  calculate 

displacements,  which  are  used  to  update  the  strain  field. 

2)  The  elastic  constants  (which  depend  on  the  current  state  of  damage)  are  used  to 

calculate  a  new  stress  field  based  on  the  updated  strain  field. 

3)  The  damage  (crack  growth)  field  is  updated  based  on  the  new  stress  field. 

4)  Return  to  step  (1)  for  another  time  increment. 

Part  (a)  of  the  damage  mechanics  is  used  in  step  (3)  to  calculate  the  increase  in  crack 
growth  associated  with  each  change  in  the  stress  field.  We  have  this  part  of  the  problem  in 
control.  Part  (b)  of  the  damage  mechanics  is  used  in  step  (2)  tc  calculate  the  stress  field 
from  the  strain  field.  This  is  the  subject  of  current  research. 


6 


The  difficulty  implementing  this  algorithm  is  that  the  elastic  constants  required  in 
step  2  are  not  simple  functions  of  the  damage,  but  depend  upon  whether  new  damage  has 
been  done  by  the  strain  increment.  The  effective  elastic  constant  is  less  when  the  cracks  are 
actively  extending  as  illustrated  by  the  stress-strain  cu^e  in  Fig.  3.  At  stresses  below  the 
fracture  initiation  stress  nonlinear  attenuation  produces  hysteresis  in  the  stress-strain  curve 
but  the  fracture  spectrum  remains  unchanged.  At  stresses  above  the  initiation  stress  the 
effective  elastic  modulus  is  lower  when  the  cracks  are  growing  than  when  they  are  simply 
sliding.  The  effective  elastic  modulus  during  subsequent  loadings  at  the  same  stress  will  be 
larger  as  illustrated  in  the  figure. 

Another  problem  in  modeling  the  effective  elasticity  of  the  damaged  rock  is  caused 
by  the  fact  that  the  crack  growth  associated  with  a  spherical  source  is  largely  radial  (parallel 
to  the  largest  principal  stress)  and  this  produces  an  axial  elastic  anisotropy  in  each  element. 
We  thus  have  to  deal  with  five  elastic  constants: 

Er  =  radial  Young's  modulus 
Et  =  transverse  Youngs  modulus 
vn  =  radial  -  transverse  Poisson's  ratio 
Vt,-  =  transverse  -  radial  Poisson's  ratio 
G  =  shear  modulus. 

The  elastic  constant.!  can  be  found  from  the  fracture  energy  relc.  se  rate  using  the  following 
expression: 


=  K2(L)  =  p2  r)C(L) 

~  E/L)  ~  2b  c)L  (8) 
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where 


G(L)  =  strain  energy  release  rate 

L  =  current  crack  length 

Ki(L)  =  stress  intensity  factor 

Er(L)  =  radial  Young's  modulus  which  we  seek 

C(L)  =  compliance 

b  =  specimen  width 

p  =  load 

Referring  to  Figure  4,  we  see  that  C=u/p,  p=sA,  u=eh,  and  Er=s/e,  so 


C  -  — h— 
FrA 


Substituting  (9)  into  (8)  gives 


K2(L)  =  -  V 

2  b  oL 


where  V=  hA  =  volume  per  crack.  Integrating  (10)  gives 


(9) 


(10) 


The  factor  2b/V  can  be  written  in  terms  of  the  crack  density  as 

2^  =  2 Nvm 

and  the  crack  density  Nv  can  be  written  in  terms  of  the  initial  damage  Dc  as 


(ID 


N  v 


=  3DoM  V3 
4  7t  laa' 


An  analytic  expression  for  Ki(L)  is  given  in  Asho/  and  Samnr.s  (1990)  which  can  be 


directly  integrated  as  in  (1 1)  to  find  Er(L). 

The  transverse  modulus,  Et ,  can  be  estimated  using  eqn.  (4)  with 
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Ki  =  a  w  tan  (r^) 

(Tada,  19X5)  for  the  transversely  loaded  array  of  cracks  shown  in  Fig  5.  Integrating  in  this 

case  gives: 

Ht  =  Ho  exp  J  2^v2  In  |  cos  (^f  )|  [  (5) 

The  shear  modulus  is  unaffected  by  crack  growth  in  this  geometry,  and  the 
Poisson's  ratios  can  he  estimated  from  the  crack  opening  displacements  (which  are  also 
functions  of  the  stress  intensity  factors). 

Figure  6  shows  the  uniaxial  stress-strain  curve  calculated  using  eqn  (11)  compared 
with  measurements  for  Berea  Sandstone  (This  work  was  done  in  collaboration  with  Randy 
Martin  and  Xiaoming  fang  at  New  Kngland  Research).  At  low  stresses,  the  negative 
curvature  of  the  data  is  caused  by  the  closing  of  cracks  and  collapse  of  pores  w'hich  are  not 
included  in  the  model.  However,  the  positive  curvature  of  the  stress-strain  curve  near 
failure  is  associated  with  the  accumulation  of  damage  and  is  well  modeled,  as  is  the  post- 
failure  decrease  in  strength.  Figure  7  shows  stress  strain  curves  for  the  different  rock  types 
studied  by  Ashby  and  Sammis  ( 19X9).  Fig.  X  shows  the  effect  of  different  levels  of  initial 
damage  on  the  elastic  behavior.  Fig.  9  explores  the  effect  of  confining  pressure,  and  Fig.  10 
shows  the  effect  ot  a  fluid  in  reducing  the  coefficient  of  friction  on  the  preexisting 
fractures.  We  are  currently  planning  further  experiments  to  test  these  predictions. 

DISCI  SSI  OX 

The  problem  of  using  the  damage  mechanics  mixlel  in  numerical  source  simulation 
cikIcn  has  been  reduced  to  the  problem  of  finding  the  stress  strain  behavior  of  a  damaged 
medium  This  is  not  a  simple  problem  because  the  effective  elastic  constant  (the  tangent 
modulus)  i>  not  a  simple  (unction  ot  crack  damage  (as  it  is.  for  example,  in  the  Biidianskv 
and  ( X'onncll  i  l!>  ;('i  iow  stiain  theory)  but  also  depends  upon  whether  an  increment  in 
sires',  produces  an  lncrcmcni  m  crack  growth.  1  he  problem  is  further  complicated  bv 

<) 


anisotropy  introduced  by  the  directional  (radial)  growth  of  fractures  in  an  explosive  stress 
field.  A  method  has  been  developed  to  calculate  the  effective  elastic  constants  from  the 
stress  intensity  factors  derived  by  Ashby  and  Sammis  (1990).  Preliminary  experimental 
tests  of  the  theory  look  promising.  Efforts  are  currently  underway  to  build  the  damage 
mechanics  into  existing  numerical  source  simulation  codes. 
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Figure  1.  Comparison  of  three  1000  pound  chemical  explosions  detonated  at  depths  of 
42,  106,  and  217  meters  in  a  limestone  quarry  (from  McEvilly  and  Johnson,  1989).  Fig. 
la  compares  the  first  arrivals  on  the  vertical  component  at  one  location.  The  pulses  have 
been  scaled  to  have  the  same  amplitude.  Fig.  lb  compares  the  amplitude  densities  of  the 
isotropic  moment  rate  tensors  for  the  three  events  while  Fig.  lc  gives  the  amplitude  density 
spectra  of  the  pulses  in  Fig.  la. 
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Figure  2.  Schematic  diagram  of  the  non-linear  zones  around  an  underground  nuclear 
explosion.  Fig.  2a  is  from  Bishop,  1963.  Fig.  2b  illustrates  the  hydrodynamic  radius  rh, 
the  damage  radius  rtj,  and  the  elastic  radius  rC|  which  are  discussed  in  the  text. 


STRESS 


igure  3.  Stress-strain  curve  illustrating  the  differences  in  behavior  above  and  below  the 
fracture  initiation  stress  (which  is  a  function  of  the  current  crack  damage).  At  the  damage 
radius  r  the  stresses  have  fallen  below  the  initiation  stress  of  the  preexisting  fractures  in  the 
emplacement  medium. 
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Figure  5.  Geometry  used  to  estimate  the  reduction  in  the  transverse  Young’s  modulus 
associated  with  crack  growth. 
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Figure  6.  Uniaxial  stress-strain  curve  for  Berea  sandstone  showing  a  comparison 
between  theory  and  experiment.  Note  that  the  reduction  in  effective  modulus  near  failure 
well  modeled. 
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Figure  7.  Comparison  of  the  uniaxial  stress-strain  behavior  of  the  various  rock  types 
studied  by  Ashby  and  Sammis  (1990).  Note  that  the  more  brittle  rocks  have  a  more 
unstable  post-failure  regime. 
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Figure  8.  Effect  of  the  initial  damage  on  the  uniaxial  stress-strain  behavior.  Note  that  the 
strength  decreases  and  the  post-failure  regime  becomes  more  stable  with  increasing  initial 
damage. 
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Figure  9.  Effect  of  the  confining  pressure  on  the  uniaxial  stress-strain  behavior.  Note 
that  the  strength  increases  and  the  post-failure  regime  becomes  more  stable  with  increasing 
confining  pressure. 
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Figure  10.  Effect  of  the  fluid- wetting  of  the  preexisting  flaws  on  the  uniaxial  stress- 
strain  behavior.  Note  that  the  strength  decreases  and  the  post-failure  regime  becomes  more 

stable  as  the  coefficient  of  the  friction  u  is  reduced  by  wetting. 
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